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We compute the total cross-section for direct Higgs boson production in hadron collisions at 
NNLO in perturbative QCD. A new technique which allows us to perform an algorithmic evaluation 
of inclusive phase-space integrals is introduced, based on the Cutkosky rules, integration by parts and 
the differential equation method for computing master integrals. Finally, we discuss the numerical 
impact of the 0{o? a ) QCD corrections to the Higgs boson production cross-section at the LHC and 
the Tevatron. 



I. INTRODUCTION 

The Higgs boson is currently the only missing particle in the minimal Standard Model (SM) of electroweak inter- 
actions. Its discovery will be one of the final steps toward the experimental verification of the SM, and will provide 
useful input for detailed studies of the mass generation mechanism and for physics beyond the SM. 

Direct searches at LEP restrict the Higgs boson mass to be greater than 114.1 GeV [1], while a global fit to 
precision electroweak measurements [2] favors a value around 90 GeV. In addition, the requirement that the SM 
remains perturbative up to relatively high energy scales sets an upper bound at approximately 1 TeV [3] . Although 
the above evidence is not completely conclusive, it indicates a relatively light Higgs boson which could be observed 
at either the Tevatron or the LHC. At both of these facilities, gluon fusion through top-quark loops is expected to be 
the dominant Higgs production mechanism. All other channels, such as vector boson fusion qq — > Hqq and associated 
Higgs production qqi — > HW, are suppressed by about an order of magnitude (see Ref. [4] for a review). We therefore 
focus upon the process gg — > H in this paper. 

The theoretical estimates of the cross-section for the Higgs boson production via gluon fusion, based on computations 
through to the next-to-leading order (NLO) in perturbative QCD, turn out to be insufficient. The leading-order (LO) 
cross-section is proportional to a^(/i 2 ), and for this reason exhibits a strong dependence on the choice of the scale \i. 
Including the 0(a s ) corrections [5,6] decreases the scale dependence, but the cross-section increases by a very large 
amount, approximately 70%. It is therefore important to evaluate the next order in the perturbative expansion, since 
this is the only way to enhance the credibility of the theoretical predictions. 

To compute the cross-section to next-to-next-to- leading order (NNLO), we must combine: the matrix elements 
for the 0(a 2 ) virtual corrections to gg — > H\ the matrix elements for the 0(a s ) virtual corrections to gg — > Hg, 
qg — > Hq, and qq — ► Hg; and finally the tree-level matrix elements for the processes gg — > Hgg, gg — > Hqq, qg — > Hqg, 
qq — + Hgg, and qq — > Hqq. For the inclusive cross-section we must integrate over the loop-momenta in the virtual 
amplitudes and the phase-space of the real particles in the final state. Both real and virtual corrections are divergent 
in four dimensions. We regularize the amplitudes using conventional dimensional regularization (d = 4 — 2e), and 
remove the ultraviolet divergences by renormalizing in the MS scheme. The remaining divergences arise from initial 
state collinear radiation and are absorbed into the parton distribution functions, yielding a finite cross-section. 

The calculation can be simplified substantially by considering the limit where the Higgs boson is much lighter 
than twice the mass of the top-quark. In this limit, the top-quark loops are replaced by point-like vertices. The 
corresponding effective Lagrangian is known to provide a satisfactory description of the cross-section for a light Higgs 
boson at NLO [5-7]. 

In the heavy top-quark limit, the NNLO contributions to the direct Higgs production cross-section are topologically 
similar to the 0{a 2 s ) corrections to the Drell-Yan process which have been calculated in the past [8]. The phase- 
space and loop integrals required for the calculation of the Higgs boson production cross-section could in principle be 
obtained in a similar fashion. However, such an approach is impractical for the Higgs production cross-section due to 
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the larger number of Feynman diagrams with a considerably more complicated tensor structure. For a problem of this 
complexity a highly automated algorithm which treats virtual and real corrections in a unified manner is desirable. 

It is well known how to construct algorithms which in principle can perform multi-loop integrations. First, one 
can employ the method of integration by parts (IBP) [9] in order to reduce the number of integrals involved in such 
computations. Algorithms which find the solutions of IBP identities in a process and topology independent manner are 
available [11,12]. After the application of IBP, a small number of remaining integrals which are not reducible further 
(master integrals), must be evaluated explicitly. Powerful techniques such as the differential equation method [10,11] 
and the Mellin-Barnes integral representation [13] can then be employed to derive an expansion of the master integrals 
in e. The above methods provide general purpose tools for the evaluation of virtual corrections. However, similar 
methods do not exist for computing phase-space integrals; they are usually calculated manually and on a case by case 
basis. In this paper we present an algorithmic procedure for evaluating phase-space integrals, based on the Cutkosky 
rules [14], integration by parts [9] and the differential equation method [11]. 

Partial results for the NNLO corrections to the Higgs boson production cross-section are available in the literature. 
The NNLO virtual corrections were computed in [15] by Harlander. The "soft" part of the cross-section at NNLO 
was derived in [16,17] by extracting contributions to the partonic cross-sections that are singular when the partonic 
center of mass energy y/l equals the mass of the Higgs boson ran- Recently, Harlander and Kilgore [18] obtained an 
excellent approximation to the complete NNLO result by expanding the phase-space integrals around the kinematic 
point s = m 2 H . In this paper we present the full analytic result for the NNLO corrections to the Higgs boson production 
cross-section. In our derivation we do not need to resort to an expansion around a special kinematic point and our 
expressions are therefore valid for an arbitrary ratio m 2 H /s. 

The paper is organized as follows. In Section II we briefly review the effective Lagrangian for describing gluon 
interactions with the Higgs boson. We also introduce our notations and present all basic formulae and definitions for 
the total cross-section. In Section III we describe our method for solving multi-particle phase-space integrals in an 
algorithmic fashion and illustrate its application with a few typical examples. We present the analytic expressions for 
the renormalized partonic cross-sections ij — > H + X in Section IV. In Section V we discuss the impact of the 0(a 2 s ) 
corrections on the Higgs boson production cross-section at the Tevatron and at the LHC. We present our conclusions 
in Section VI. Some useful formulae, including the complete list of master integrals, are collected in the Appendix. 



II. EFFECTIVE LAGRANGIAN 



The Higgs boson interaction with gluons is a loop induced process and is therefore sensitive to all colored particles 
which get their masses through the Higgs mechanism. In this paper we restrict ourselves to the Standard Model where 
the top-quark contribution dominates. 

Although the Born cross-section is known as a function of the top mass m t and the Higgs boson mass run, it is 
much harder to obtain the exact analytic dependence of the cross-section on the mass of the top-quark in higher 
orders of perturbation theory. However, since it is most probable that the Higgs boson is light, it is sufficient to work 
in the infinite top-quark mass limit. 

For Higgs boson masses in the range 100 — 200 GeV, we can describe the Higgs gluon interaction by introducing 
the effective Lagrangian [5-7] 



(1) 



where G" is the gluon strength tensor, H is the Higgs field and v w 246 GeV is the Higgs boson vacuum expectation 
value. The Wilson coefficient C\, defined in the MS scheme, is [19] 
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where a s {y) is the MS strong coupling constant, nf — 5 is the number of active flavors and L t — \og([i 2 /m 2 ). 

It is expected that the effective Lagrangian of Eq. (1) is a valid approximation to the Higgs gluon interaction for 
small values (mjj < 2m t ) of the Higgs boson mass. It can be checked that at leading order and for mn ~ 150 GeV, 
the effective Lagrangian approximation is accurate within 5%, whereas for m H <~ 200 GeV, the accuracy drops to 
10%. The precision of the approximation improves for the Higgs boson production cross-section computed at NLO 
accuracy [6] . The effective Lagrangian description therefore seems accurate in the entire range of phenomenologically 
interesting Higgs boson masses, and we adopt it for the calculation of the NNLO corrections. 
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The effective Lagrangian approach separates the short-distance (<~ m^ 1 ) and the long-distance (~ Tn^ ) scales, 
simplifying the calculation of the Higgs boson production cross-section. For example, the original one-loop triangle 
diagram of the gluon-gluon Higgs interaction vertex at LO is now replaced by the simple tree-level vertex derived 
from Eq.(l). The effective Lagrangian approach also yields the correct H ggg and Hgggg interaction vertices, as a 
consequence of gauge invariance. In addition, in the limit of vanishing fermion masses, there is no direct interaction 
between massless quarks and the Higgs boson. 

The partonic cross-sections for the production of the Higgs boson, up to NNLO in perturbation theory, receive 
the following contributions: a) virtual corrections to gg — > H, up to 0(a 2 ); b) virtual corrections to single real 
emission processes gg — > Hg 1 qg — > Hq, qg — > Hq, qq — > Hg, up to 0(a s ); and c) double real emission processes 
gg — ► Hgg, gg — > Hqq, qg — > ifqg, — ► i?g<7, — > 99 — > to LO. The effective Lagrangian of Eq.(l) and 

the corresponding matrix elements should be renormalizcd in the MS scheme by a global renormalization factor [19]: 
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where /?o an d /?i are the two first coefficients of the QCD /3-function: 
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This additional renormalization, together with the standard renormalization of the strong coupling constant removes 
all ultraviolet divergences in the cross-section. However, since we work in the approximation of massless colored 
partons in the initial state, the total cross-section is not finite even after the ultraviolet renormalization has been 
performed. The remaining singularities are associated with the collinear radiation off the colliding partons. It is well 
known that these singularities factorize and can be removed by renormalizing the parton distribution functions in a 
manner consistent with the DGLAP evolution equation. 

The general factorization formula for the cross-section of the Higgs boson production from the collision of two 
hadrons /ii(pi) and h 2 {p 2 ), is 



cr/n+ha— ff+x = ^ / dx 1 dx 2 fj; hl \x 1 )f ( j h2 \x 2 )a i3 ^ H (m 2 H ,x 1 x 2 s), 

ij o 
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where f\ h \x) is the standard distribution function for a parton i in the hadron h, Uij is the partonic cross-section for 
i + j — > H + X and s = (p\ +P2) 2 is the square of the total center of mass energy of the hadron hadron collision. Using 
dimensional analysis we can write the partonic cross-sections in terms of the single dimensionless variable z = m 2 H /s, 

(Tij - l/v 2 g{m 2 H /s) (6) 

with s = sxix 2 - Introducing x = m 2 H /s, we then rewrite Eq. (5) in the form 

<T hl+ H^H+X =*E [fi Hl) ® ^ ® K-WA)] (7) 
ij 

where the standard convolution ® is defined as 

1 

[/i®/2](i) = y dx 1 dx 2 fi(x 1 )f 2 {x 2 )6(x - xix 2 ). (8) 


The collinear singularities are factored out from the partonic cross-sections with the following procedure. Denoting 
the unrenormalized (in the sense of collinear singularities) partonic cross-sections by cry , the MS renormalizcd partonic 
cross-sections &ij are implicitly given by: 

[<Jij(x)/x] = Yy°ki{z)/z] O T kt ® F y , (9) 



where the kernels F^ arc 
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The standard space-like splitting functions P,j [20,21] are listed in the Appendix. We can easily solve Eq.(9) for 
&ij(x)/x = Pij{x) order by order in a s . It is convenient to introduce a matrix notation and rewrite Eq. (9) as 



p = T 1 (gipOT, 



(11) 



where p is the matrix of partonic cross-sections in flavor space and T is the matrix with components (x) as in 
Eq.(10). We then write 



T = S(l-x)U-^T 1+ (^) T 2 + 0(a 3 s ), 



(12) 

where the matrix U has the components JZy = Si g Sj g and I\2 can be read off from Eq. (10). Inverting Eq.(ll) to 
obtain 



P = [r^^gpotrp 1 
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and expanding p in a s 



we find 
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Having derived the finite partonic cross-sections <7y, we must convolute them with the MS parton distribution 
functions fi to obtain the total hadronic cross-section: 



° hl+h ^H + X = [jf ® ® (*«(*)/*)] (*). 



(15) 



We present our results for the partonic cross-sections &ij(z) in Section IV. In Section V we use Eq.(15) to calculate 
the Higgs boson production cross-section at the Tevatron and at the LHC. 



III. METHOD 



In this Section we describe the method employed to compute the partonic cross-sections to NNLO. At this order 
we must calculate three distinct contributions: 

• double- virtual: the interference of the Born and the two-loop amplitude as well as the self-interference of the 
one-loop amplitude for gg — ► H, 




+ 148 terms; 



• real- virtual: the interference of the one-loop and the Born amplitudes for gg — > Hg, gq — > Hq, and gq — > Hq, 




+ 635 terms; 



• double-real: the self-interference of the Born amplitudes for gg — > Hgg, gg — > Hqq, gq — > Hgq 7 gq — > Hgq, 
qq —>■ Hqq, and qq — > -ff<7<7, 
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+ 594 terms. 



The above interference terms are produced in a form convenient for further evaluation using the QGRAF package [22] 
for generating Feynman graphs. 

In the following subsections we briefly describe the available techniques for evaluating virtual corrections and explain 
our method for integrating over the phase-space of the final state particles. 



A. Virtual corrections 



There currently exists a general method which permits the systematic evaluation of multi-loop virtual corrections. In 
order to calculate a multi-loop amplitude we must first reduce the number of Feynman integrals. The hypergeometric 
structure of Feynman integrals guarantees that simple algebraic relations between various scalar integrals exist, making 
such a reduction possible. One method of producing these relations is the integration by parts (IBP) technique [9]. 
In cases where the system of equations is not complete, it can be supplemented with additional identities that exploit 
the Lorentz invariance (LI) of scalar integrals [11]. 

In general, IBP and Lorentz Invariance (LI) identities relate integrals of differing complexity. For example, it 
is possible that a single IBP equation relates an integral with an irreducible scalar product to integrals with no 
irreducible scalar products or to integrals with fewer propagators. A typical situation however, involves multiple IBP 
and LI identities relating several equally complicated integrals to a set of simpler ones. In such cases, every integral 
must be written exclusively in terms of simpler ones and, eventually, expressed in terms of a few "master" integrals 
which cannot be reduced further. Unfortunately, finding recursive solutions of the IBP and LI identities is tedious, 
and may be impossible in complicated cases. Also, a separate treatment of each different topology in a Feynman 
amplitude is required. Consequently, the whole procedure becomes increasingly cumbersome with the introduction of 
more kinematic variables and loops. 

We may alternatively consider a sufficiently large system of explicit IBP and LI equations which contains all the 
integrals that contribute to the multi-loop amplitude of interest. It should then be possible to solve the system 
of equations in terms of the master integrals [11,12] using standard linear algebra elimination algorithms. In this 
approach, the number of loops, the topological details, and the number of kinematic variables, affect only the size of 
the system of equations and the number of terms in each of the equations; they have no bearing on the construction 
of the elimination algorithm. This in principle allows us to express any multi-loop amplitude in terms of master 
integrals. 

One possible elimination algorithm has been proposed by Laporta [12]. This algorithm exploits the fact that 
Feynman integrals can be ordered by their complexity; for example they can be arranged according to the number 
of irreducible scalar products and the total number and powers of propagators. This observation distinguishes the 
IBP and LI systems of equations from algebraic systems with no intrinsic ordering, and it becomes possible to solve 
them iteratively, starting with the simpler equations and progressing to more complicated. We use a variant of this 
algorithm, implemented in FORM [23] and MAPLE [24] . 

After the reduction we must compute the analytic expansion in e of the master integrals. The coefficients of the 
expansion are typically expressed in terms of polylogarithms whose rank and complexity depends on the number of 
loops and kinematic variables of the integral in question. The Mcllin-Barnes representation [13] and the differential 
equation method [11] can be used to evaluate master integrals explicitly. 



B. Reduction of phase-space integrals 



In this subsection we extend the application of the above techniques to calculate phase-space integrals for inclusive 
cross-sections. To the best of our knowledge the method we present is new, however a somewhat related discussion 
has been given earlier in [25]. 

To illustrate our method, we consider the following double-real contribution at NNLO: 
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Using the Cutkosky rules [14], we can replace the delta-functions in the above integral by differences of two propagators: 

1 1 



2inS (p — m 



The r.h.s. of Eq. (16) is now equal to a forward scattering diagram: 
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where a cut propagator should be replaced by the r.h.s. of Eq.(17). 

We have exchanged the square of a Born amplitude for a two-loop diagram, in contrast to the usual application of 
the Cutkosky rules. We do this in order to utilize IBP and LI relations between multi-loop integrals. The phase-space 
integrals can then be evaluated in the same algorithmic fashion as the multi-loop integrals. 

We begin our calculation by summing over the colors and spins of the external particles in the cut two-loop integral 
on the right hand side of Eq.(18). The original diagram is then expressed in terms of a large number of scalar 
two-loop integrals to which the same cutting rules apply. Crucially, we can use the IBP method to reduce the cut 
scalar integrals. This is a consequence of the fact that the delta-function in Eq. (17) is represented in a very simple 
manner by the difference of two propagators with opposite prescriptions for their imaginary parts. We derive the 
IBP equations by integrating over total derivatives which act on the propagators of the cut scalar integrals. The 
prescription for the imaginary part of the two propagators in the r.h.s. of Eq. (17) is irrelevant for the differentiation. 
Therefore the IBP relations for the two descendants of these two terms have the same form as the IBP relations for 
the original integral without the cut. It is then allowed to commute the application of IBP reduction algorithms with 
the application of the Cutkosky rules. 

After the IBP reduction, the original phase-space integral is expressed in terms of a small number of master integrals 
cut through the same three propagators as the initial diagram 1 : 



= A l 



(19) 



During the reduction, integrals with one or more of the cut propagators eliminated are produced. From Eq. (17) 
we observe that such terms do not contribute to the original phase-space integrals. Therefore, we can immediately 
discard them simplifying the reduction process. 

A similar procedure can be applied to the virtual-real contributions. In this case, since we perform the phase-space 
integration over two final state particles, the resulting master integrals should be cut through two of the propagators: 




Bo 




(20) 



In order to have a unified algorithm for all three types of interferences, we treat the double-virtual corrections as 
integrals with a single cut through the propagator of the Higgs boson. 



1 Bold lines represent a massive Higgs propagator. Normal lines denote massless scalar propagators. 
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Finally we must evaluate the master integrals as a series expansion in e = (4 — d)/2. Since each of the cut master 
integrals represents a well-defined phase-space integral, we could compute them using brute force techniques similar 
to the ones described in Ref. [8]. However, we can instead utilize the IBP reduction algorithm in order to produce 
a set of coupled first order differential equations [11] that the master integrals satisfy. It is simpler to solve the 
differential equations than to reinstate the delta-functions for the cut propagators and perform the integrations over 
the phase-space. 



C. Evaluation of phase-space master integrals 

To explain how the system of differential equations for the master integrals is obtained, let us consider a two-loop 
scalar integral with a single Higgs boson propagator 

2 f d d k d d l 1 

X V> m «)-J (2^(2^ \V- m \\ v ...A^- (22) 

By differentiating with respect to m 2 H we obtain: 

dl(s,m 2 H ) f _d^^H 1 

dm% V J (2nY (2ny [k 2 _ m y+i A *. . . . A ^ ' { 6 > 

After applying the IBP algorithm, we can rewrite the r.h.s of the last equation in terms of the master integrals {A^}, 
yielding 
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By identifying X = Xi, we derive a closed system of differential equations for the master integrals, 
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These differential equations can be solved up to a constant in terms of logarithms and generalized Nielsen polyloga- 
rithms, order by order in e. This constant is obtained by evaluating the master integral at a specific kinematic point. 
This is typically simpler and can often be avoided using general arguments, as shown in an example below. 

As an explicit example we discuss the calculation of the master integrals for the real-virtual corrections. The IBP 
reduction produces six master integrals which depend on two variables: the mass of the Higgs boson, ran, and the 
square of the sum of the incoming momenta, s = (pi + P2) 2 ■ Note that the integrals depend on a single Mandelstam 
variable, since they correspond to forward scattering diagrams with the same incoming and outgoing momenta. It is 
convenient to express the master integrals in terms of the dimensionless ratio z = m 2 H /s. We can further simplify our 
results by setting s = 1. The full dependence on s can be restored by simple dimensional analysis. 

Three of these master integrals are combinations of the one-loop massless bubble integral and the two-body phase- 
space integral and can be easily evaluated, yielding 
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have a more complicated dependence on z and we compute them using the method of differential equations. As 
an example, we discuss the differential equation for the first integral in Eq.(30): 





This differential equation is of the form 
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and has the general solution 

f(z) = J' dxa(x) (c + jT dxp{x)e- f dx ' a ^ . 
Using Eq.(33) and the expressions for the two boundary master integrals in Eqs. (27,28), we derive 
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We compute the value of the integral at a specific kinematic point in order to determine the constant C. A convenient 
choice is the threshold for Higgs production, z = 1, where the integral vanishes: 




(2 = 1) 



(35) 



The value of the integral at z — 1 can be inferred without an explicit calculation by observing that, in the z — > f limit 
the two particle phase-space scales like (1 — z), and the one-loop triangle diagram on the r.h.s. of the cut scales like 
(f — z)°. Requiring that Eq. (34) vanishes at z = 1, one finds: 
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We now evaluate the integrals in Eq.(34). First, by changing the integration variables, we isolate the singularity at 
x = 0: 
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The integrand on the r.h.s. can then be expanded in e, yielding generalized Nielsen poly logarithms, 
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Substituting this result in Eq.(34) and truncating the series at the order where polylogarithms of rank n + p > 3 start 
to appear, we arrive at the result: 
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We repeat the same procedure for the differential equations for the two remaining integrals in Eq.(30), where the 
master integral we have just calculated enters as a boundary term. It is important to extract the singular behavior 
of the master integrals around z = 1 before expanding in e. This is essential since terms of the form (1 — z)~ 1+ac in 
the cross-section are expanded in e in terms of "plus" distributions, 
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facilitating the cancellation between real and virtual soft and collinear singularities prior to integration over z. 

Finally, we apply the same technique to compute the double real master integrals. Explicit formulae for the required 
master integrals are given in the Appendix. 



IV. PARTONIC CROSS-SECTIONS 

In this section we present analytic expressions for the partonic cross-sections i + j — > H + X of Eq. (15). We write 
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and a s is the MS strong coupling constant evaluated at the scale /i r = ma- For simplicity, the factorization scale is 
also set equal to the mass of the Higgs boson /i/ = ma- 







At leading order we find 
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At next-to-leading order there arc contributions from the gluon-gluon, quark-gluon and quark-antiquark channels: 
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The main result of this paper is the next-to-next-to-leading order corrections, which we separate according to their 
dependence on the number of quark flavors: 



For the gluon-gluon channel we find 
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2 1 + s Ll3( "" } _ 2 1^ L l2 (x)ln(x) 



9 (2a; 4 - 15a; 2 - 10a; - 5) 
2 1 + x 

27(3x 2 + 2x + 1) 



Li 2 (-x) ln(x) 



9 (59 + 177a; 2 - 116x 3 + 59x 4 - 118a;) 



1-x 



ln(x) ln 2 (l-x) 



1 + x 



T . / u \ 9(6-llx 3 + 18x 2 -12x + 6x 4 ), 2 , 
Li 2 (-x) ln(l + x) + ln(x) ln(l - x) 



1-x 



+ 2 



9(3-8x 3 + 3x 4 -6x + 9x 2 ) 



+ 



1-x 

9 (8x 4 + 16x 3 + 33x 2 + 22x + 11) 
2 1 + x 

9 (4x 4 + 8x 3 + 27x 2 + 18x + 9) 
4 1 + x 



Li 2 (x) ln(l x) — 



3 (7x - 7x 3 + 4 + 18x 2 - 17x 4 + 9x 5 ) 



1-x 2 



ln 3 (x) 



C2 l n (i + x ) _ m{x2 + X + l) \ \ 2 ( x ) ln(l + x) 



1 + x 



ln(l + x) ln%x) + (-21 + ^-x 2 - 18x + 



63 

T 



33 



5 ) ln(l + x) ln(x) 



2687 2027, , . 

x + — — )ln(l-x) 



27 (3x 2 + 2x + 1) , , . 3 (-280x 3 + 143x 4 + 394x- 289 + 21x 2 ) T . , . 
+ y _ In (l + x ) ln(x)- i _ L l2 (x) 

, 63 2 ,„ 33 ,. T . . . , 2559 , 1079 
+ (-21 + yx 2 - 18x + _x 3 )Li 2 (-x) + ( — x 3 + —^-s 

3 (374x 4 - 389x+ 154 + 699x 2 - 827x 3 ) , 2 . . . „ , ojo , oo nml , 2 , 
--^ — '- ln 2 (x) + (330x 3 - 348x 2 + 381x - 297) ln 2 (l - x) 

3 (-1180x 3 + 641 - 1238x+ 1227x 2 + 605x 4 ) , , . , , 2x , 

+ -^ '- ln(x) ln(l - x) - 72(2 - x + x 2 )x ln 3 (l - x) 

1 (4318x 4 - 6955x 3 + 6447x 2 - 561 lx + 2333) , . . 3 (495x 4 - 886x 3 + 564x 2 - 200x + 16) , 
"8 —x ln(x) + 4 T— x C2 

9(6x + 18x 2 + 2 + 10x 5 - 6x 3 - 19x 4 ) A , , . 9 (-48x 3 + 23x 4 - 46x + 3 + 69x 2 ) A , . 
+ — - 2 (2 ln(x) - _ ( 2 ln(l - x) 



+ 2 



9 (-36 - 15x 4 - 52x + 19x 2 + 13x 3 + 33x 5 ) 



1-x 2 



7539 



24107 
^8" 



-x^ + 



22879 
IT* 



18157 



(47) 



(48) 
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and 



. (2)F . 1189 5 A 5 A 2 T . . 10 

A ^ =(-144 + 6 C3 -3 C2+ 3 W -^ 3 



ln(l - a;) 



+2 



+ 



"ln(l - 


x) 2 ~ 




1 - 


a; 


+ 


47 2 


25 


12* 4 


1 

~6* + 


1 


1 


6 + 


12 



,'31 1 65 2 , . 
h ' ~6 X+ 6 + U X ) Sl2 W + 



Li 2 (x) ln(x) 



1 - X 

31 
12 
1 
6 

11 1 2 

12 12" ' 24 T 



1-x 



1 17 



x Li 3 (x) 



1 2 1 

12* + 6^ 



x 2 Li 2 (x) ln(l -»)+—- — x + — x 2 ln(l - x) ln(x) In 



6 6 

C 2 ln(l - x) - 4x(l + x)( 2 ln(x) 
x)ln^ 



V x 



5 m 3 Z' 17 2 7 !\ ^ / 34 . 2 2 8 16\ 2/1 . A . 

--x(l + x)ln 3 x+ (^-_ a; 2 -- a; --JC 3 + g-x 3 + -x 2 --x + -J (ln 2 (l-x)-C 2 ) 

2 (21x 2 + 7x + 25x 4 + 17 - 61x 3 ) 



1-x 



. 785 3 83 2 49 461\ , , 
ln(x) ln(l - x) + ( — x 3 - -x 2 + -x - — ) ln(l - x) 



1 (-351x 3 + 117x 2 + 68 + 132x 4 + 52x) , 2 . . 1 (227x 3 + 68 + 4x 4 - 302x + 21x 2 ) T . ,„ 
+ ™ ; In ( x )+^ i Li 2 (l-x) 



+ 



72 1-x 
1 (333x 2 + 2384x 4 - 598x - 3041x 3 + 1282) 



36 



216 1-x 
For the quark-gluon channel we obtain 



ln(x) 



8887^ 
648 



1267 2 497 
-x — — — x 



432 



216 



12923 
1296 ' 



(49) 



170 



338 119 



x 2 Li 3 (x) + (4x + 4 + 2x 2 )Li 3 (-x) + (16 + 8x 2 + 16x)S 12 {-x) 



367 367 



367 



614 269 2 7 4\ „ , s , „ o , , Nr< , o, fou, uu! 2 ™i \, a,, 

; 1 ~~ x ~ ~ x ~ yj Sl2{x) + { ~ 2x ~ 4 " 4x ^ x ) + (-27 + -54* - ) ln ( J - ' ; 



+ (2 + x 2 - 2x) ln(l - x) - 



446 



214 281 
~ + ~9~ 



x 2 ln(x) - (8 + 4x 2 + 8x) ln(l + x) Li 2 (x 



/o o . o,, / (1 + a;) \ T . . . / 115 2 230 230 \ , , N1 2 , 
+ (8 + 8x + 4x 2 ) ln ( V - ; J Li 2 (-x) + ( — x 2 + — x - — J ln(x) ln 2 (l - x) 

107 107 2 107 V 2, x ( 145 2 71 „V , 

T + "IT ""J ln ( * } ln(1 _ x) + r^T" _ T7 X ~ 2 J ln(x) 

+ (-3x 2 - 6 - 6x) ln(l + x) ln(x) 2 + (4x + 4 + 2x 2 ) ln(l + x) 2 ln(x) 

4 , 74 11 2 166\ , . /2605 146 74 , 79 ,\ 
+ I -y 7 * 3 - lf X ~~9 ~^2n) U ^ + {-^-— X+ 27 X ~6 X ) L12(X) 



^. x + ^ + Sx 3 -72] ln 2 (l-x) + 
18 12 J 

,'113 3 244 13 2 31 V 2 . , 
M^ 3 + -x- y x 2 - y ln 2 (x) + 



121 

l8~ ; 



326 



4 3 74 

7^X 6 - —X 

27 9 



826 5935\ . „ 
- —x + J ln(x) ln(l - x) 

^x 2 -^)ln(l + x)ln(x) 



, 59 2 H8 118V . , /140 128 2 52 \ , 

+Cz ( -~9 X + — x -—) ln(1 ~ x) +C2 + y ' ln(a 



„ , n „ N , 392 3 49 2 23671 
+C 2 (12 + 12x + 6x 2 ) ln(l + x) + ( -— x 3 - -x 2 + 



106xj ln(l 



1985 

ToV 



800 



12209 616 ,V . . / 292 3 82 16 2 221 \ , 

1 x 3 ln(x) + x 3 x H x 2 H C2 

162 81 ) ' \ 27 3 3 27 / s 



92 



+ (-18x+10x^ + -lC 3 - m4 



210115 1537 



+ 



486 



■x- 5 + ■ 



16465 
~L62 



-x + 



2393 
648 S 



(50) 



and 



/I 7 1 1\, 2 „ x ( 38 19 2 29 \ , , N 209 265 
A « = ( 18* " 9* + 9 ln (1 - X) + -27* + 27^ + 27 ^ " ^ + 162 
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+ (j- - + -* - -x 2 J ln(x) - + -x -j HI x) + ^x 2 + (^ 2 - + g j In (*). 

For the scattering of two identical quarks we obtain: 

Am4 /368 104 2 400\ . . 32, , 2ri , . 4 .„ 2 „ 2 , 
™ = ["27 * + "27~ + "27~J 3(X) " ¥ (X + } 12(aj) " 27 ( + } (X) ( ~ 

A 1 f-1 

- — (a; + 2) 2 ln 3 (x) - —(19 + 5.x 2 + 17x)Li 2 (x) ln(x) - — (x + 3)(1 - x) ln 2 (l - x) 
27 27 9 

+ ^(x + 3)(1 - x) ln(x) ln(l - x) + ^-(26x - 18 + 9x 2 ) ln 2 (x) - ^(-6 + x 2 + 4x)Li 2 (x) 
o Zi ( y 

4,, „„w., N1 x /8. 118 248 46 2 \ , . , 

+ -(5x + 17)(1 - x) ln(l - x) + (^-(x + 2) 2 C 2 - — + ~^x + -x 2 J ln(x) 

/ 8 o 16 16 \ , (Vo 32 8 ,V 4 ,„„ 
+ - ^ + 27* Us + y - y x - -x 2 C2 - ^(27x + I60)(l - x), 



and 



Al 2)F = 0. 



</</ 



For the scattering of distinct quarks we find 

A T = f ( x + 2 ) 2 ( Li 3W - 5i2(a:)) - l(x + 2) 2 ln(x)Li 2 (x) - ±(x + 2) 2 ln 3 (x) 

8 32 16 

--(4a; - 6 + x 2 )Li 2 (x) - — (x + 3)(1 - x) ln 2 (l - x) + — {x + 3)(1 - x) ln(x) ln(l - x) 

9 9 3 

+ ^(x 2 + 4x - 3) ln 2 (x) + ^( 2 (x + 2) 2 ln(x) + Ubx + 17)(1 - x) ln(l - x) 
y y o 

2 / 16 32 8 \ 2 

+ -(29x 2 + 44x - 59) ln(x) + ( — - —x - -x 2 J C2 - g (Hx + 105)(1 - x), 



A ( r = 0. 

99' 



and 

,(2)F 
99' 

Finally, for the quark-antiquark channel the NNLO contribution is 

^-(-^-§*K*K-^-P-i') &( -*> 

32 32 4 4 

+ — (x + 2) 2 Li 3 (x) - — (x + 2) 2 5i 2 (x) - — (x + 2) 2 ln 3 (x) + -(2 + 2x + x 2 ) ln(l + x) ln 2 (x) 
■ ' y z, ( y 

~(2 + 2x + x 2 ) ln 2 (l + x) - I (a; + 2) 2 Li 2 (x) + |(2 + 2x + x 2 )Li 2 (-x)^ ln(x) 

-i|(2 + 2x + x 2 )Li 2 (-x) ln(l + x) + |^(1 - x)(13x 2 - 35x - 14) ln 2 (l - x) 

27 81 

16 8 
- — (1 - x)(37x 2 - lOlx - 44) ln(x) ln(l - x) - — (44x 3 + 39x - 81x 2 + 27) ln 2 (x) 

81 81 

16 8 
+— x(x + 6x 2 + 2) ln(l + x) ln(x) + — (42x - 87x 2 + 12 + 10x 3 )Li 2 (x) 
27 81 

+— x(x + 6x 2 + 2)Li 2 (-x) - — (1 - x)(384x 2 - 967x - 75) ln(l - x) + ( -— x 2 - — - — x 
27 81 V 27 27 27 

+ (§„ + 2 ft 2 + ^ - ^ - - ,»(,) - A<2 + 2, + ,n,l + ,) 

and 
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Ag>' = m ( i -*) + (- 



81 



64 
"2T 



64 3 16 80 , 
81* 27 + 27 X ) lnW 



243 



(1-x)(41.t 2 -88.t + 23). 



(57) 



The above results are valid if the renormalization and factorization scales are equal to the mass of the Higgs boson, 
(i r = (if = run- It is easy to restore the complete functional dependence of the partonic cross-sections on these scales 
using the fact that the total hadronic cross-section is independent of them. 

We first find the dependence of the partonic cross-sections on a scale (i which is equal to both the renormalization 
and factorization scales (i r = (if = (i. To do so, we restore the full (i dependence in Eq.(15): 



(58) 



Since the physical cross-section <Jh 1 +h 2 ^H+x does not depend on (i, 

o d 



and the derivatives of the structure functions fl with respect to (i can be determined from the DGLAP evolution 
equation: 



0. 



(59) 



we derive the relation 



ijk 



(hi) 



a s ((i) 



Pi 



ik 



(S^Hx^) = ^ [Pj ® /,(z,/i)] (x), 



(a kj (z,(i)/z) + M 2 ^2 (^v( z ^)/ z ) + (vik(z,(i)/z) ^Atl p kj 



A(i 2 



p(h 2 ) 



(60) 



(61) 



This equation should hold for arbitrary (i and x; therefore, the expression in the square brackets should be identically 
zero for any choice of i and j, yielding the following "evolution equation" for the cross-section: 



d a s {n) 



Pik <® {crkj {Z, (l)/z) + (&ik (z, (i)/z) ® Pkj 



(62) 



We can solve Eq.(62) order by order in a s ((i) using the partonic cross-sections at (i = inn as the boundary condition. 
Having obtained the explicit dependence of the partonic cross-sections on (i, 



(63) 



we can find their dependence on two independent renormalization and factorization scales by expressing the strong 
coupling constant a s ((i = (if) through a s {(i r ). 

We have checked that our expressions for the partonic cross-sections, derived by explicitly evaluating the Feynman 
amplitudes with their full scale dependence, are in agreement with Eq. (62). Our results are also in complete agreement 
with Ref. [18] where the first sixteen terms of an expansion of the partonic cross-sections in 1 — x were computed. 
In the limit, x — > 0, only the leading logarithmic corrections are known [26]. We can easily reproduce this result by 
expanding our formulae for partonic cross-sections around x = 0. 



V. NUMERICAL RESULTS 



We can now discuss the numerical impact of the NNLO corrections on the Higgs boson production cross-section 
at the LHC and the Tevatron. To calculate the cross-section we must convolute the hard scattering partonic cross- 
sections of Section IV with the appropriate parton distribution functions. For a self-consistent calculation at NNLO, 
we need the parton distribution functions at a given factorization scale at the same order. At present, the NNLO 
evolution of the distribution functions can not be performed since the required three-loop splitting functions arc not 
known. Nevertheless, a significant number of moments of the splitting functions is available [27], and this information 
can be combined with the known behavior at small x [28] , to obtain a useful approximation for the NNLO splitting 
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functions [29]. In Ref. [30] this approach was used to determine the NNLO MRST parton distribution functions; we 
use these for the numerical evaluation of the Higgs boson production cross-section. 

To demonstrate the convergence properties of the perturbative series for the hadronic cross-section, we present 
the LO, NLO and NNLO results for both the LHC and the Tevatron. In order to improve on the heavy top-quark 
approximation, we normalize the cross-section to the exact Born cross-section with the full m t dependence. We use 
the mode = 1 parton distribution functions (see Ref. [30] for the notation). For the NNLO set, this mode provides 
the "average" of two extreme cases, the so-called fast and slow evolutions. For the evaluation of the strong coupling 
constant we use LO, NLO and NNLO running accordingly, with the Z-po\e values used in the parton distribution 
functions as initial conditions (see [30] for details). The total cross-section for the LHC is shown in Fig.l. We note that 
the NNLO cross-section does not vary significantly if we choose a different mode for the MRST parton distribution 
functions; the observed changes are less than 1%. 



a{pp^H + X) [pb] , y/l = 14 TeV 

100 r 1 . 1 




1 1 1 ' ' 1 

100 150 200 250 300 



tuh, GeV 

FIG. 1. The Higgs boson production cross-section at the LHC at leading (dotted), next-to- leading (dashed-dotted) and 
next-to-next-to-leading (solid) order. The two curves for each case correspond to fx r = Hf = m-H /2 (upper) and fj, r = fj,f = 2m,H 
(lower). 

From Fig.l we observe that the scale dependence of the Higgs production cross-section at NNLO in the range 
/i r = fif = m#/2 — 2m,H is approximately 15%; this is a factor of two smaller than the NLO scale dependence and a 
factor of four smaller than the LO variation. Despite the scale stabilization, the corrections are rather large; the NLO 
corrections increase the LO cross-section by about 70%, and the NNLO corrections increase it further by approximately 
30%. The K factor, defined as the ratio of the NNLO cross section and the LO cross-section at fi r = fif = tuh, is 
approximately two. In Fig. 2 we plot the values of the Higgs production cross-section at the Tevatron. The NNLO K 
factor is approximately three, and the residual scale dependence is approximately 23%. These numerical results for 
the NNLO Higgs boson production cross-section are in excellent agreement with the corresponding results in [18]. 



a(pp - 


>H + X) 


pb], V^=2TeV 



























100 150 200 250 300 

run, GeV 

FIG. 2. The Higgs boson production cross-section at the Tevatron at leading (dotted), next-to- leading (dashed-dotted) and 
next-to-next-to-leading (solid) order. The two curves for each case correspond to fx r = Hf = mij/2 (upper) and fj, r = /// = 2m,H 
(lower). 
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A few remarks concerning the magnitude of the corrections are appropriate. Despite the fact that the mass of 
the Higgs boson is much smaller than the total center of mass energy, the production cross-section is dominated by 
partonic processes with s ~ m%, This is because the gluon gluon luminosity is a rapidly decreasing function of the 
partonic center of mass energy. The agreement between our numerical results based on the complete expressions of 
Section IV with the approximate results of [18], where an expansion in 1 — m 2 H /s was employed, demonstrates this 
indirectly. 

The dominance of the threshold region renders resummation methods applicable [31,32]. However, threshold dom- 
inance should also affect the cross-section estimates based on fixed order calculations where there is freedom in the 
choice of the factorization scale. Since the production process is dominated by the region x — > 1, the appropriate 
factorization scale should be parametrically smaller than the mass of the Higgs boson; choosing a factorization scale 
near the Higgs boson mass may not capture the essential physics of the process. 

We illustrate this point by considering the NLO correction to the Higgs production cross-section. Concentrating 
on the gluon-gluon subprocess, and keeping the most singular terms in the x — * 1 limit, we can write 

^ )(x) = (v){(y + 6C2 ) ^-*)- 6 

It is obvious from the above expression that if the dominant contribution to the integrated cross-section comes from 
the region x ~ 1, then choosing /i = ra# leaves large logarithmic corrections of the form log(l — x) in the hard 
scattering cross-section. To avoid this problem, we should choose /i <~ mij(l — x), which is parametrically smaller 
than the mass of the Higgs boson. While it is not possible to use an x-dependent factorization scale without resorting 
to a full resummation program, in the fixed order calculation we can attempt to do this on average. This choice 
decreases the NNLO corrections and the Higgs boson production cross-section increases as compared to conventional 
choice of the scales, \i T = (if = tuh- 

a(pp -> H + X) [pb], s/s = 14 TcV 

100 ,— 
90 - 
80 
70 . 
60 ... 
50 - 
40 .... 
30 
20 
10 - 

20 30 40 50 60 70 80 90 100 50 100 150 200 250 

/x, GeV /x, GeV 

FIG. 3. The Higgs boson production cross-section at the LHC at leading (dotted), next-to- leading (dashed-dotted) and 
next-to-next-to-leading (solid) order as the function of factorization and renormalization scale a. The mass of the Higgs boson 
is 115 GeV for the left and 275 GeV for the right plot. 

We demonstrate this behavior with two examples in Fig. 3, where we plot the production cross-sections for mu = 
115 GeV and m# = 275 GeV. We equate the renormalization and factorization scales and vary the factorization 
scale from \x = 15 GeV up to the mass of the Higgs boson. These plots illustrate that for smaller values of /j, the 
NLO cross-section increases more rapidly than the NNLO cross-section, and the difference between the NLO and the 
NNLO results becomes smaller. Therefore, the convergence of the perturbative series is improved for smaller values 
of the factorization scale. 

If we adopt this argument and restrict our analysis to small p, we find a Higgs production cross-section of 55 ± 5 pb 
for rriH = 115 GeV, a somewhat larger value than obtained with the conventional scale choice /x = tuh • It is interesting 
that recent studies [32] of the threshold resummed cross-section for Higgs boson production, matched to the NNLO 
calculation, detect a similar increase as compared to fixed order calculations with /x = tuh- 



l-x 



•In 



m 2 H {l — x)' 



+ 



(64) 
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VI. CONCLUSIONS 



In this paper, we have studied the Higgs boson production cross-section in hadron-hadron collisions. The main 
contribution to the hadronic cross-section originates from gluon-gluon fusion, which we have computed at NNLO 
(O (a*)) in perturbative QCD. The other partonic production channels, qg^>H + X,qq^>H + X,qq^H + X 
and qq' — > H + X, were also studied to order a A s . 

We have presented explicit analytic expressions for the partonic cross-sections valid within the heavy top-quark 
approximation. Finally, we have calculated the cross-section for direct Higgs boson production at the Tevatron and 
the LHC by performing a numerical convolution of the partonic cross-sections with the MRST 2002 NNLO set [30] 
of parton distribution functions. The residual scale dependence of the NNLO cross-section is approximately <~ 15% 
for the LHC and ~ 23% for the Tevatron. The NNLO if-factors are fairly large for both the LHC and the Tevatron. 
Nevertheless, the cross-section increases less from NLO to NNLO than from LO to NLO, indicating a slow convergence 
of the perturbative expansion. 

While this calculation was in progress, Harlandcr and Kilgore [18] obtained an approximation of the partonic cross- 
sections by expanding around the Higgs boson production threshold. Our results for the partonic cross-sections are 
in complete agreement with their expansion. Also, we find a very good agreement of our numerical results for the 
NNLO Higgs production cross section with the results reported in [18]. 

In Section V, we have argued that it is more appropriate to choose smaller values of the factorization scale /i 
than the conventional choice /i = ran- Then, the NNLO corrections decrease, indicating a better convergence of the 
perturbative series. Moreover, with this choice, the fixed order results are in better agreement with recent estimates 
of the cross-section based on threshold resummation [32] . 

We have suggested a method for the algorithmic evaluation of inclusive phase-space integrals. This method combines 
the Cutkosky rules with integration-by-parts reduction algorithms to achieve a systematic reduction of phase-space 
integrals to a few master integrals. We have also shown how to compute these master integrals using differential 
equations produced with the IBP reduction algorithms. 

The techniques discussed in this paper can be used to compute higher order corrections to other inclusive processes of 
direct phenomenological interest. We are also confident that this approach can be generalized to enable the calculation 
of differential distributions. In fact, a connection between phase-space integrals with a modified measure, such as 
the integrals that appear in the evaluation of invariant mass, energy and angular distributions, and loop integrals 
with unconventional propagators exists. This connection can be used to automate the calculation of differential 
distributions following the lines of Section III. The above ideas will be the subject of more detailed studies in future 
work. 

Acknowledgments Wc would like to thank Robert Harlander and Bill Kilgore for comparisons with their results 
of Ref. [18]. We are grateful to Frank Petriello for careful reading of the manuscript and invaluable suggestions. This 
research was supported by the DOE under grant number DE-AC03-76SF00515. 
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APPENDIX 



In this Appendix, some formulae used in our calculation are given. 



A. Splitting functions 

We first summarize the formulae for the space-like splitting functions. 
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B. Master Integrals 

Below we list the master integrals required in this calculation. We denote s = p\ 2 = (Pi + P2) 2 , z — m 2 H /s and for 
simplicity we set s = 1. 

1. Double-Virtual 

The master integrals for the double-virtual corrections can be expressed in terms of Gamma functions with the 
exception of the cross-triangle which has been calculated in [33]. For completeness, we list their e-expansion below: 
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2. Real-Virtual 



For the real- virtual contributions we find six master integrals: 
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3. Double-Real 



We find the following master integrals for the double-real contributions: 
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